Exercise solutions: Section 8.11

Ex 1

Figure 8.31 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers.

1. Explain the differences among these figures. Do they all indicate that the data are white noise?

• The figures show different critical values (blue dashed lines).

• All figures indicate that the data are white noise.

2. Why are the critical values at different distances from the mean of zero? Why are the autocorrelations different in each figure when they each refer to white noise?

• The critical values are at different distances from zero because the data sets have different number of observations. The more observations in a data set, the less noise appears in the correlation estimates (spikes). Therefore the critical values for bigger data sets can be smaller in order to check if the data is not white noise.

Ex 2

A classic example of a non-stationary series is the daily closing IBM stock price series (data set ibmclose). Use R to plot the daily closing prices for IBM stock and the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.

ggtsdisplay(ibmclose)

The time plot shows the series “wandering around”, which is a typical indication of non-stationarity. Differencing the series should remove this feature.

ACF does not drop quickly to zero, moreover the value \(r_1\) is large and positive (almost 1 in this case). All these are signs of a non-stationary time series. Therefore it should be differenced to obtain a stationary series.

PACF value \(r_1\) is almost 1. All other values \(r_i, i>1\) are smaller than the critical value. These are signs of a non-stationary process that should be differenced in order to obtain a stationary series.

Ex 3

For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.

1. usnetelec

autoplot(usnetelec)

There is no need for a Box-Cox transformation in this case.

usnetelec %>% diff() %>% autoplot()

2. enplanements

autoplot(enplanements)

enplanements %>% BoxCox(lambda=0) %>% diff(lag=12) %>% diff() %>% autoplot()

3. visitors

autoplot(visitors)

visitors %>% BoxCox(lambda=0) %>% diff(lag=12) %>% diff() %>% autoplot()

Ex 4

For the enplanements data, write down the differences you chose above using backshift operator notation.

A seasonal (lag 12) difference, followed by a first difference: \[ (1-B) (1-B^{12})\]

Ex 5

myts <- readxl::read_excel("data/retail.xlsx", skip=1)[,"A3349873A"] %>%
  ts(frequency=12, start=c(1982,4))

First we will try logs as before, followed by a seasonal difference.

myts %>% log() %>% diff(lag=12) %>% ggtsdisplay()

That still looks a little non-stationary, so I will take another difference at lag 1.

myts %>% log() %>% diff(lag=12) %>% diff() %>% ggtsdisplay()

That looks stationary now. So we used a log, a seasonal difference, and a first difference. Other time series from retail data set may require a different set of transformation/differencing choices.

Ex 6

Use R to simulate and plot some data from simple ARIMA models. a. Use the following R code to generate data from an AR(1) model with \(\phi_{1} = 0.6\) and \(\sigma^2=1\). The process starts with \(y_1=0\).

ar1 <- function(phi, n=100)
{
  y <- ts(numeric(n))
  e <- rnorm(n)
  for(i in 2:n)
    y[i] <- phi*y[i-1] + e[i]
  return(y)
}

2. Produce a time plot for the series. How does the plot change as you change \(\phi_1\)?

Some examples of changing \(\phi_1\)

autoplot(ar1(0.6))

autoplot(ar1(0.95))

autoplot(ar1(0.05))

autoplot(ar1(-0.65))

3. Write your own code to generate data from an MA(1) model with \(\theta_{1} = 0.6\) and \(\sigma^2=1\).

ma1 <- function(theta, n=100)
{
  y <- ts(numeric(100))
  e <- rnorm(100)
  for(i in 2:100)
    y[i] <- theta*e[i-1] + e[i]
  return(y)
}

4. Produce a time plot for the series. How does the plot change as you change \(\theta_1\)?

autoplot(ma1(0.6))

autoplot(ma1(0.95))

autoplot(ma1(0.05))

autoplot(ma1(-0.8))

5. Generate data from an ARMA(1,1) model with \(\phi_{1} = 0.6\), \(\theta_{1} = 0.6\) and \(\sigma^2=1\).

arma11 <- function(phi, theta, n=100)
{
  y <- ts(numeric(100))
  e <- rnorm(100)
  for(i in 2:100)
    y[i] <- phi*y[i-1] + theta*e[i-1] + e[i]
  return(y)
}
autoplot(arma11(0.6,0.6))

6. Generate data from an AR(2) model with \(\phi_{1} =-0.8\), \(\phi_{2} = 0.3\) and \(\sigma^2=1\). (Note that these parameters will give a non-stationary series.)

ar2 <- function(phi1, phi2, n=100)
{
  y <- ts(numeric(100))
  e <- rnorm(100)
  for(i in 3:100)
    y[i] <- phi1*y[i-1] + phi2*y[i-2] + e[i]
  return(y)
}
autoplot(ar2(-0.8,0.3))

7. Graph the latter two series and compare them.

See graphs above. The non-stationarity of the AR(2) process has led to increasing oscillations

Ex 7

Consider wmurders, the number of women murdered each year (per 100,000 standard population) in the United States.

1. By studying appropriate graphs of the series in R, find an appropriate ARIMA(\(p,d,q\)) model for these data.

wmurders %>% ggtsdisplay()

The graphs suggest differencing of the data before applying ARMA models.

wmurders %>% diff() %>% ggtsdisplay()

The ACF and PACF graphs suggest that the resulting time series might be an AR(2) or MA(2) due to the significant spike at lag 2 in both graphs. So either an ARIMA(0,1,2) or an ARIMA(2,1,0) appear appropriate. As we need to make a choice, I’ll use an ARIMA(0,1,2) in what follows.

2. Should you include a constant in the model? Explain.

A constant would imply a drift in the original data which does not look correct, so we omit the constant.

3. Write this model in terms of the backshift operator.

ARIMA(0,1,2): \((1-B)y_t = (1 + \theta_1B + \theta_2B^2) \varepsilon_{t}\).

4. Fit the model using R and examine the residuals. Is the model satisfactory?

(fit <- Arima(wmurders, c(0,1,2)))

## Series: wmurders 
## ARIMA(0,1,2) 
## 
## Coefficients:
##          ma1    ma2
##       -0.066  0.371
## s.e.   0.126  0.164
## 
## sigma^2 estimated as 0.0422:  log likelihood=9.71
## AIC=-13.43   AICc=-12.95   BIC=-7.46

checkresiduals(fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,2)
## Q* = 9.77, df = 8, p-value = 0.28
## 
## Model df: 2.   Total lags used: 10

• No obvious patterns are visible on the graph of the residuals
• Ljung-Box test suggests that the residuals are white noise
• The residuals look close enough to normally distributed

The model is satisfactory.

5. Forecast three times ahead. Check your forecasts by hand to make sure that you know how they have been calculated.

(fcast <- forecast(fit, h = 3))

##      Point Forecast  Lo 80  Hi 80  Lo 95  Hi 95
## 2005         2.4585 2.1952 2.7217 2.0558 2.8611
## 2006         2.4771 2.1169 2.8373 1.9262 3.0280
## 2007         2.4771 1.9793 2.9749 1.7157 3.2385

To check the forecasts by hand we use the model written in terms of the backshift operator: \[ (1-B)y_t = (1 + \theta_1B + \theta_2B^2) \varepsilon_{t}, \] where \(y_t\) is the observed time series and the coefficients are \(\theta_1 = -0.066\) and \(\theta_2 = 0.371\). It can be rewritten as: \[ y_t = y_{t-1} + \varepsilon_t -0.066\varepsilon_{t-1} + 0.371\varepsilon_{t-2}. \] The last observation is 2.5894 and the last few residuals are -0.39266, -0.34376, 0.05024. So the next forecast is: \[ \hat{y}_{T+1|T} = 2.5894 + 0 -0.066 (0.0502) + 0.371 (-0.3438) = 2.4585. \] Similarly, \[ \hat{y}_{T+2|T} = 2.4585 + 0 + 0 + 0.371 (0.0502) = 2.4771 \] and \[ \hat{y}_{T+3|T} = 2.4771 + 0 + 0 + 0 = 2.4771. \]

6. Create a plot of the series with forecasts and prediction intervals for the next three periods shown.

autoplot(fcast)

7. Does auto.arima give the same model you have chosen? If not, which model do you think is better?

fit.auto <- auto.arima(wmurders)
fcast.auto <- forecast(fit.auto, h = 3)
autoplot(fcast.auto)

auto.arima picked a more complex model which assumes that the trend in the data will continue in the future. The extra level of differencing is hard to justify from the graphs, and it is unlikely that murder rates will fall indefinitely. So I’d probably prefer only one difference. We can force auto.arima to do that:

(fit.auto = auto.arima(wmurders, d=1, stepwise=FALSE))

## Series: wmurders 
## ARIMA(0,1,2) 
## 
## Coefficients:
##          ma1    ma2
##       -0.066  0.371
## s.e.   0.126  0.164
## 
## sigma^2 estimated as 0.0422:  log likelihood=9.71
## AIC=-13.43   AICc=-12.95   BIC=-7.46

Now auto.arima picks the same model that I did.

Ex 8

Consider austa, the total international visitors to Australia (in millions) for the period 1980-2015.

1. Use auto.arima() to find an appropriate ARIMA model. What model was selected. Check that the residuals look like white noise. Plot forecasts for the next 10 periods.

(fit <- auto.arima(austa))

## Series: austa 
## ARIMA(0,1,1) with drift 
## 
## Coefficients:
##         ma1  drift
##       0.301  0.173
## s.e.  0.165  0.039
## 
## sigma^2 estimated as 0.0338:  log likelihood=10.62
## AIC=-15.24   AICc=-14.46   BIC=-10.57

An ARIMA(0,1,1) with drift model is selected .

checkresiduals(fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1) with drift
## Q* = 2.3, df = 5, p-value = 0.81
## 
## Model df: 2.   Total lags used: 7

The residuals look like white noise.

fit %>% forecast(h=10) %>% autoplot()

2. Plot forecasts from an ARIMA(0,1,1) model with no drift and compare these to part a. Remove the MA term and plot again.

austa %>% Arima(order=c(0,1,1), include.constant=FALSE) %>%
  forecast() %>% autoplot()

austa %>% Arima(order=c(0,1,0), include.constant=FALSE) %>%
  forecast() %>% autoplot()

Removing the drift term makes a big difference as the forecasts no longer have a linear trend. Removing the MA(1) term makes almost no difference.

3. Plot forecasts from an ARIMA(2,1,3) model with drift. Remove the constant and see what happens.

austa %>% Arima(order=c(2,1,3), include.constant=TRUE) %>%
  forecast() %>% autoplot()

austa %>% Arima(order=c(2,1,3), include.constant=FALSE) %>%
  forecast() %>% autoplot()

Error in stats::arima(x = x, order = order, seasonal = seasonal, include.mean = include.mean,  :
  non-stationary AR part from CSS

The model fitting causes an error because the drift is really needed. What happens is that the AR terms try to handle the drift but end up going outside the stationarity region.

4. Plot forecasts from an ARIMA(0,0,1) model with a constant. Remove the MA term and plot again.

austa %>% Arima(order=c(0,0,1), include.constant=TRUE) %>%
  forecast() %>% autoplot()

austa %>% Arima(order=c(0,0,1), include.constant=FALSE) %>%
  forecast() %>% autoplot()

A stationary model with a constant has long term forecasts equal to the mean (long term here meaning after the first horizon). A stationary model without a constant has long term forecasts equal to zero.

5. Plot forecasts from an ARIMA(0,2,1) model with no constant.

austa %>% Arima(order=c(0,2,1), include.constant=FALSE) %>%
  forecast() %>% autoplot()

The second level of differencing induces a linear trend in the forecasts, but the prediction intervals become quite wide.

Ex 9

For the usgdp series: a. if necessary, find a suitable Box-Cox transformation for the data;

autoplot(usgdp)

I don’t think a Box-Cox transformation is required.

2. fit a suitable ARIMA model to the transformed data using auto.arima();

(fit <- auto.arima(usgdp))

## Series: usgdp 
## ARIMA(2,2,2) 
## 
## Coefficients:
##          ar1    ar2     ma1     ma2
##       -0.123  0.311  -0.584  -0.367
## s.e.   0.287  0.087   0.300   0.286
## 
## sigma^2 estimated as 1604:  log likelihood=-1199.6
## AIC=2409.1   AICc=2409.4   BIC=2426.4

3. try some other plausible models by experimenting with the orders chosen;

The second order differencing will induce a trend in the forecasts, which is required here, so I will look at changing only \(p\) and \(q\).

fit020 <- Arima(usgdp, order=c(0,2,0))
fit021 <- Arima(usgdp, order=c(0,2,1))
fit022 <- Arima(usgdp, order=c(0,2,2))
fit023 <- Arima(usgdp, order=c(0,2,3))
fit120 <- Arima(usgdp, order=c(1,2,0))
fit121 <- Arima(usgdp, order=c(1,2,1))
fit122 <- Arima(usgdp, order=c(1,2,2))
fit123 <- Arima(usgdp, order=c(1,2,3))
fit220 <- Arima(usgdp, order=c(2,2,0))
fit221 <- Arima(usgdp, order=c(2,2,1))
fit222 <- Arima(usgdp, order=c(2,2,2))
fit223 <- Arima(usgdp, order=c(2,2,3))
fit320 <- Arima(usgdp, order=c(3,2,0))
fit321 <- Arima(usgdp, order=c(3,2,1))
fit322 <- Arima(usgdp, order=c(3,2,2))
fit323 <- Arima(usgdp, order=c(3,2,3))

4. choose what you think is the best model and check the residual diagnostics;

The best according to the AICc values is the ARIMA(2,2,1) model.

checkresiduals(fit221)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,2,1)
## Q* = 9.36, df = 5, p-value = 0.095
## 
## Model df: 3.   Total lags used: 8

The residuals pass the Ljung-Box test, but the histogram looks like it has heavier tails than Gaussian.

5. produce forecasts of your fitted model. Do the forecasts look reasonable?

fit221 %>% forecast(h=20) %>% autoplot

These look reasonable.

6. compare the results with what you would obtain using ets() (with no transformation).

usgdp %>% ets %>% forecast(h=20) %>% autoplot

The ETS point forecasts are higher than the ARIMA forecasts, and the ETS forecast intervals are wider.

Ex 10

Consider austourists, the quarterly number of international tourists to Australia for the period 1999–2010.

1. Describe the time plot.

autoplot(austourists)

lambda = BoxCox.lambda(austourists)

BoxCox.lambda function suggests using Box-Cox transformation with parameter \(\lambda=0.064\). This value is very close to \(0\) and moreover Box-Cox transformation with parameter \(0\) has another meaning: it is log transformation which ensures that the forecasted values transformed back to the original scale (using exponent transformation) will be positive. Since number of tourists visiting Australia cannot be negative it makes sense. Therefore we will just use logs of the original data.

austourists %>% log() %>% ggtsdisplay()

The ACF and PACF graphs suggest seasonal differencing of the data before applying ARMA models.

tsdisplay(diff(log(austourists), lag=4))

The ACF and PACF graphs suggest that the resulting time series can be ARMA(0,1)(0,1) or ARMA(0,1)(1,0) or ARMA(1,0)(0,1) or ARMA(1,0)(1,0) models. Therefore the suggested ARIMA model for log(austourists) data are ARIMA(0,0,1)(0,1,1) or ARIMA(0,0,1)(1,1,0) or ARIMA(1,0,0)(0,1,1) or ARIMA(1,0,0)(1,1,0). A constant should be included in the model since the differenced time series is raised significantly above zero. The constant will lead to a linear trend in the forecasts. Such behaviour is strongly supported by the data.

fit1 = Arima(austourists, c(0,0,1), c(0,1,1), include.drift = TRUE, lambda = 0)
fit2 = Arima(austourists, c(0,0,1), c(1,1,0), include.drift = TRUE, lambda = 0)
fit3 = Arima(austourists, c(1,0,0), c(0,1,1), include.drift = TRUE, lambda = 0)
fit4 = Arima(austourists, c(1,0,0), c(1,1,0), include.drift = TRUE, lambda = 0)

The best model out of these four can be chosen using AICc:

c(fit1$aicc, fit2$aicc, fit3$aicc, fit4$aicc)

## [1] -167.84 -156.36 -170.40 -159.63

The lowest AICc is for model ARIMA(1,0,0)(0,1,1) with drift, therefore it is the preferred model for forecasting.

auto.arima(austourists, lambda = 0)

## Series: austourists 
## ARIMA(1,0,0)(2,1,0)[4] with drift 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##         ar1    sar1    sar2  drift
##       0.433  -0.710  -0.237  0.013
## s.e.  0.118   0.126   0.134  0.002
## 
## sigma^2 estimated as 0.00415:  log likelihood=85.65
## AIC=-161.3   AICc=-160.27   BIC=-150.51

auto.arima finds the same model according to AICc.

The best model written using backshift operator and without it

Using backshift operator the model can be written as: \((1 - \phi_1 B) (1 - B^4) y_t = c + (1 + \theta_1 B^4) \epsilon_t\)

Without the backshift operator the model can be written as: \(y_t - \phi_1 y_{t-1} - y_{t-4} + \phi_1 y_{t-5} = c + \epsilon_t + \theta_1 \epsilon_{t-4}\) or: \(y_t = \phi_1 y_{t-1} + y_{t-4} - \phi_1 y_{t-5} + c + \epsilon_t + \theta_1 \epsilon_{t-4}\)

The last form allows forecasting value \(y_t\) assuming that values \(y_{t-1}\), \(y_{t-4}\), \(y_{t-5}\) and \(\epsilon_{t-4}\) are known. Value \(\epsilon_t\), although, should be estimated as zero.

2. What can you learn from the ACF graph?

3. What can you learn from the PACF graph?

4. Produce plots of the seasonally differenced data \((1 - B^{4})Y_{t}\). What model do these graphs suggest?

5. Does auto.arima() give the same model that you chose? If not, which model do you think is better?

6. Write the model in terms of the backshift operator, then without using the backshift operator.

Ex 11

Consider usmelec, the total net generation of electricity (in billion kilowatt hours) by the U.S. electric industry (monthly for the period January 1973 – June 2013). In general there are two peaks per year: in mid-summer and mid-winter.

1. Examine the 12-month moving average of this series to see what kind of trend is involved.

2. Do the data need transforming? If so, find a suitable transformation.

3. Are the data stationary? If not, find an appropriate differencing which yields stationary data.

4. Identify a couple of ARIMA models that might be useful in describing the time series. Which of your models is the best according to their AIC values?

5. Estimate the parameters of your best model and do diagnostic testing on the residuals. Do the residuals resemble white noise? If not, try to find another ARIMA model which fits better.

6. Forecast the next 15 years of electricity generation by the U.S. electric industry. Get the latest figures from the EIA to check the accuracy of your forecasts.

7. Eventually, the prediction intervals are so wide that the forecasts are not particularly useful. How many years of forecasts do you think are sufficiently accurate to be usable?

Data set usmelec

usmelec

##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 1973 160.22 143.54 148.16 139.59 147.40 161.24 173.73 177.37 156.88 154.20
## 1974 157.56 142.75 150.34 142.31 153.81 156.44 178.25 174.12 152.47 152.20
## 1975 164.62 147.35 155.76 146.50 153.53 162.72 177.06 179.93 155.44 155.19
## 1976 178.61 156.97 164.47 153.47 157.66 173.67 186.69 186.64 165.24 164.01
## 1977 196.66 162.95 169.44 157.12 169.60 181.03 199.17 196.36 176.50 166.65
## 1978 198.11 173.75 173.46 160.01 175.55 188.59 202.95 206.66 185.80 176.01
## 1979 209.99 186.59 183.15 170.26 178.41 186.98 202.52 205.10 180.97 179.95
## 1980 200.30 188.96 187.75 169.02 176.07 189.75 217.06 215.63 191.70 178.76
## 1981 206.76 179.86 185.83 172.84 178.14 203.02 220.66 210.64 187.05 181.56
## 1982 209.69 180.55 187.97 172.88 177.48 186.45 210.87 205.89 180.88 173.17
## 1983 195.87 172.72 182.77 170.67 174.72 191.37 220.45 230.19 195.82 183.14
## 1984 216.92 189.81 200.39 181.38 192.55 209.97 221.53 229.53 195.41 191.14
## 1985 228.15 198.49 195.25 185.17 197.12 205.68 227.00 226.29 202.71 195.00
## 1986 217.76 192.58 197.12 186.37 197.65 215.33 242.95 225.40 206.91 197.96
## 1987 223.04 194.28 202.13 189.79 206.41 225.91 248.20 247.88 213.22 203.22
## 1988 238.19 217.18 214.29 196.30 208.70 233.07 257.74 267.93 220.39 210.81
## 1989 246.77 233.78 241.95 222.87 234.82 250.92 273.46 275.42 242.79 235.67
## 1990 255.19 229.50 244.76 229.77 241.77 268.99 287.45 290.65 258.36 244.37
## 1991 269.21 228.29 240.56 227.80 254.87 269.10 294.65 292.16 256.03 245.38
## 1992 267.77 239.51 247.73 233.41 242.41 261.08 293.62 281.93 259.93 244.99
## 1993 271.02 248.01 261.25 234.69 244.33 275.36 312.23 311.45 264.03 250.55
## 1994 289.77 249.17 258.00 240.64 252.75 294.16 311.26 307.61 266.26 256.53
## 1995 279.77 252.31 261.34 244.74 264.29 286.26 330.42 345.78 277.57 263.98
## 1996 296.92 270.69 275.02 251.61 282.27 302.72 327.71 329.29 283.15 271.45
## 1997 300.57 258.13 272.26 258.28 272.91 299.09 344.52 332.90 301.06 284.97
## 1998 295.26 260.59 286.88 261.23 299.64 328.90 361.94 357.37 318.92 284.45
## 1999 315.81 274.82 298.14 280.72 300.10 328.92 376.54 364.03 305.52 283.94
## 2000 327.99 294.17 301.58 285.58 322.95 339.05 356.53 368.67 312.45 289.45
## 2001 332.49 282.94 300.71 278.08 300.49 327.69 357.61 370.53 306.93 294.73
## 2002 319.94 281.83 302.55 289.85 307.68 341.02 381.54 374.59 331.28 307.06
## 2003 341.99 299.25 304.32 285.76 307.55 328.69 374.40 381.82 323.14 306.74
## 2004 346.55 314.28 308.81 290.56 327.38 345.08 377.33 368.44 335.62 312.45
## 2005 343.12 298.50 317.46 289.56 315.06 363.67 402.27 404.94 350.22 316.40
## 2006 328.66 307.33 318.73 297.86 330.62 364.26 410.42 407.76 332.06 321.57
## 2007 353.53 323.23 320.47 303.13 330.20 362.75 393.23 421.80 355.39 332.62
## 2008 363.00 325.11 324.63 305.87 325.25 373.11 402.90 388.99 338.06 318.55
## 2009 354.99 300.89 310.60 289.54 311.31 347.66 372.54 381.22 327.40 307.04
## 2010 360.96 319.74 312.17 287.80 327.94 375.76 409.73 408.88 346.05 307.92
## 2011 363.11 313.29 318.71 302.40 323.63 367.73 418.69 406.54 337.96 308.73
## 2012 340.92 310.15 309.04 295.94 337.53 361.51 416.51 396.11 334.74 312.16
## 2013 348.64 309.60 325.37 298.26 322.12 356.40                            
##         Nov    Dec
## 1973 148.14 153.60
## 1974 150.07 160.01
## 1975 153.03 169.63
## 1976 169.35 184.14
## 1977 167.39 184.59
## 1978 176.39 192.10
## 1979 177.77 188.97
## 1980 178.77 195.85
## 1981 175.79 195.83
## 1982 173.60 184.96
## 1983 183.17 212.56
## 1984 190.60 200.23
## 1985 192.65 219.49
## 1986 196.65 213.79
## 1987 200.48 220.74
## 1988 209.81 232.99
## 1989 234.25 274.45
## 1990 231.24 255.77
## 1991 241.66 254.09
## 1992 244.16 267.34
## 1993 252.10 272.17
## 1994 251.87 269.52
## 1995 261.48 285.55
## 1996 269.12 284.25
## 1997 271.12 296.36
## 1998 267.14 297.99
## 1999 271.10 295.17
## 2000 284.67 319.01
## 2001 278.93 305.50
## 2002 296.29 324.83
## 2003 297.87 331.68
## 2004 302.10 341.95
## 2005 306.12 348.10
## 2006 309.16 336.28
## 2007 314.10 346.29
## 2008 310.05 343.90
## 2009 296.63 350.51
## 2010 306.01 362.12
## 2011 304.12 335.75
## 2012 305.55 334.33
## 2013

plot(usmelec); lines(ma(usmelec, 12), col = "red", lwd = 2)

The data show strong up-trend which might be considered as linear in long run.

The seasonal pattern which size increases together with the level of the trend strongly suggests using Box-Cox transformation before forecasting.

lambda = BoxCox.lambda(usmelec)
plot(BoxCox(usmelec, lambda))

BoxCox.lambda function suggests using Box-Cox transformation with parameter \(\lambda=-0.57\).

The data is not stationary (it contains trend), therefore differencing is required to make it stationary for ARMA forecasting.

tsdisplay(diff(BoxCox(usmelec, lambda), lag = 12))

ACF and PACF graphs suggest that the data is stationary. Suggested models are: ARIMA(1,0,0)(0,1,1) ARIMA(2,0,0)(0,1,1) ARIMA(3,0,0)(0,1,1)

fit1 = Arima(usmelec, c(1,0,0), c(0,1,1), lambda = lambda)
fit2 = Arima(usmelec, c(2,0,0), c(0,1,1), lambda = lambda)
fit3 = Arima(usmelec, c(3,0,0), c(0,1,1), lambda = lambda)

c(fit1$aicc, fit2$aicc, fit3$aicc)

## [1] -5002.6 -5032.7 -5058.1

The lowest AICc is for model ARIMA(3,0,0)(0,1,1).

Residuals analysis

tsdisplay(fit3$residuals)

Box.test(fit3$residuals, lag = 5, type = "Ljung-Box", fitdf = 4)

## 
##  Box-Ljung test
## 
## data:  fit3$residuals
## X-squared = 11.2, df = 1, p-value = 8e-04

ACF and PACF graphs as well as Ljung-Box test suggest that the residuals are not white noise. Looking at the residuals we can suggest spikes at lags 3 and 4 at ACF graph can be removed by increasing MA value in ARIMA model:

fit4 = Arima(usmelec, c(3,0,1), c(0,1,1), lambda = lambda)
fit5 = Arima(usmelec, c(3,0,2), c(0,1,1), lambda = lambda)
fit6 = Arima(usmelec, c(3,0,3), c(0,1,1), lambda = lambda)

c(fit4$aicc, fit5$aicc, fit6$aicc)

## [1] -5089.8 -5087.0 -5085.8

It appears that currently model ARIMA(3,0,1)(0,1,1) is with the lowest AICc.

The residuals for this model are:

tsdisplay(fit4$residuals)

Which are good except for the spikes at ACF and PACF graphs at lag 15. We will try to rectify it by increasing AR and MA values for seasonal and non-seasonal ARIMA parts:

fit7 = Arima(usmelec, c(3,0,1), c(1,1,1), lambda = lambda)
fit8 = Arima(usmelec, c(4,0,1), c(0,1,1), lambda = lambda)
fit9 = Arima(usmelec, c(4,0,1), c(1,1,1), lambda = lambda)

c(fit7$aicc, fit8$aicc, fit9$aicc)

## [1] -5089.1 -5087.8 -5087.2

All received AICc values are greater than for the best found model ARIMA(3,0,1)(0,1,1). Therefore ARIMA(3,0,1)(0,1,1) should be used for forecasting. The residuals although appear to be correlated for lag 15:

Box.test(fit4$residuals, lag = 15, type = "Ljung-Box", fitdf = 4)

## 
##  Box-Ljung test
## 
## data:  fit4$residuals
## X-squared = 27.8, df = 11, p-value = 0.0034

It means that the forecasting intervals, calculated using this model can be incorrect.

Forecasting 15 years ahead and comparing with existing data (from http://www.eia.gov/totalenergy/data/monthly/#electricity)

# New data taken from http://www.eia.gov/totalenergy/data/monthly/#electricity website
newData = c(306009.629,362119.081,362871.893,313126.607,318709.941,302400.724,323628.238,
         367727.015,418692.755,406511.315,337931.318,308698.504,304102.155,335740.463,
         339528.259,309389.42,309090.569,295228.216,336517.631,360826.151,414639.671,
         395699.566,334584.877,311650.774,305975.035,334635.089,348489.737,309435.119,
         325301.486,298073.696,321833.617,356224.262,393799.21,383968.062,340293.27,
         314682.729,313751.509,352356.534,377019.193,323662.011,331595.192,296766.045,
         323730.672,357419.408,384838.868,383494.071,338975.899,313971.928,317176.051)/1000
newTs = ts(newData, start = c(2010, 11), frequency = 12)

fcast = forecast(fit4, h = 15*12)
plot(fcast)
lines(newTs, col = "red", lwd = 1.5)

accuracy(fcast, newTs)

##                    ME    RMSE    MAE       MPE   MAPE    MASE      ACF1
## Training set -0.13594  7.2298 5.3642 -0.014499 2.0398 0.59571 -0.030083
## Test set      1.05573 11.0216 8.5795  0.410227 2.4113 0.95278  0.593354
##              Theil's U
## Training set        NA
## Test set       0.37767

The forecasts are good. Probably quite a few years of forecasts can be correct. Although the forecasting intervals are expanding very quickly and therefore the forecasts become less and less useful when the forecasting horizon increases.

Ex 12

For the mcopper data:

1. if necessary, find a suitable Box-Cox transformation for the data;

2. fit a suitable ARIMA model to the transformed data using auto.arima();

3. try some other plausible models by experimenting with the orders chosen;

4. choose what you think is the best model and check the residual diagnostics;

5. produce forecasts of your fitted model. Do the forecasts look reasonable?

6. compare the results with what you would obtain using ets() (with no transformation).

Data set mcopper

mcopper

##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 1960  255.2  259.7  249.3  258.0  244.3  246.8  250.6  241.3  231.0  218.7
## 1961  216.6  220.1  221.7  225.6  238.6  232.8  226.1  227.2  225.8  225.0
## 1962  226.8  231.3  231.1  230.6  230.5  230.4  230.4  230.4  230.4  230.6
## 1963  230.6  230.6  230.6  230.6  230.6  230.6  230.6  230.6  230.6  230.6
## 1964  234.2  247.8  266.1  307.7  295.8  288.8  305.4  356.9  415.0  485.1
## 1965  356.7  419.2  440.6  480.5  490.8  466.1  404.0  431.6  473.5  500.1
## 1966  599.0  668.7  668.7  679.9  592.8  604.8  559.5  426.4  402.4  454.9
## 1967  443.7  435.4  391.8  354.9  369.5  362.3  355.9  372.7  378.3  406.3
## 1968  586.4  715.8  708.0  522.4  456.5  473.0  439.9  439.9  462.3  449.7
## 1969  523.8  535.0  532.4  578.1  579.6  617.0  605.4  669.0  657.6  647.4
## 1970  677.2  690.1  730.4  725.2  665.9  607.0  567.8  527.5  519.3  475.8
## 1971  420.7  424.9  476.5  521.2  464.1  447.2  464.4  451.0  427.5  417.7
## 1972  418.9  426.9  442.2  433.5  423.3  412.4  423.2  427.0  434.2  428.7
## 1973  474.7  511.8  610.2  638.8  613.0  678.4  795.4  843.5  799.6  849.6
## 1974  913.2 1006.5 1172.0 1267.7 1190.6 1020.1  802.8  767.9  630.1  599.2
## 1975  512.4  528.7  554.5  560.5  539.7  522.5  559.3  603.7  580.2  573.1
## 1976  587.8  601.6  683.7  817.9  836.0  878.0  922.1  862.2  844.9  784.8
## 1977  814.8  833.4  881.0  830.4  797.9  763.1  724.4  665.8  685.6  682.9
## 1978  651.2  627.0  657.6  694.4  716.0  725.0  705.4  734.6  735.8  750.1
## 1979  827.1  969.7 1005.5 1012.3  935.3  888.8  802.4  883.2  953.6  967.6
## 1980 1147.9 1220.2 1045.5  937.1  888.9  858.5  916.7  877.9  857.6  846.2
## 1981  777.4  784.9  814.3  837.0  834.0  861.0  897.3  981.4  942.2  904.9
## 1982  853.6  863.8  836.3  858.5  843.6  740.4  829.9  841.0  832.7  861.4
## 1983  997.3 1074.3 1071.9 1090.0 1122.8 1098.6 1115.6 1091.0 1041.0  958.4
## 1984  978.0  991.5 1030.9 1078.1 1022.3  991.2 1007.9 1018.2 1029.3 1043.5
## 1985 1205.5 1269.8 1234.7 1212.7 1225.3 1118.0 1067.5 1025.7 1001.4  973.8
## 1986  995.4  982.8  984.4  957.3  932.3  937.0  891.7  876.5  916.1  922.9
## 1987  893.8  902.7  920.1  909.3  911.9  964.6 1052.8 1097.5 1100.7 1182.7
## 1988 1476.8 1324.3 1286.3 1216.5 1307.1 1428.8 1297.5 1296.2 1445.7 1689.4
## 1989 1912.6 1765.0 1904.1 1832.3 1679.1 1638.8 1539.1 1731.5 1834.9 1801.8
## 1990 1432.7 1390.9 1615.7 1639.6 1633.7 1510.5 1529.6 1554.4 1611.6 1409.6
## 1991 1265.1 1246.5 1322.5 1412.3 1336.8 1344.8 1353.8 1325.6 1346.1 1371.2
## 1992 1182.2 1240.5 1291.7 1260.9 1224.6 1239.1 1313.8 1297.2 1307.0 1360.3
## 1993 1472.3 1536.7 1472.3 1261.9 1159.0 1228.5 1287.5 1305.0 1219.7 1108.7
## 1994 1208.4 1261.0 1284.0 1268.1 1430.5 1549.9 1589.5 1560.0 1601.0 1586.2
## 1995 1910.5 1830.3 1826.8 1805.7 1745.9 1876.8 1927.4 1936.6 1870.3 1782.7
## 1996 1708.7 1650.2 1676.5 1713.5 1757.3 1407.7 1276.6 1294.9 1244.3 1236.2
## 1997 1469.2 1480.3 1506.4 1466.8 1540.0 1588.0 1466.7 1403.6 1315.4 1256.4
## 1998 1032.1 1014.4 1051.4 1078.1 1057.8 1005.6 1004.2  991.8  979.1  935.6
## 1999  866.7  866.5  849.6  910.2  935.6  891.3 1041.2 1025.5 1077.3 1040.1
## 2000 1124.2 1124.7 1100.6 1059.8 1183.3 1161.6 1192.3 1245.1 1365.3 1308.1
## 2001 1209.8 1214.9 1202.8 1159.3 1178.9 1147.5 1078.5 1019.0  974.4  948.5
## 2002 1049.9 1097.4 1127.9 1101.4 1092.9 1110.0 1022.2  962.4  950.0  952.5
## 2003 1018.8 1046.7 1047.4 1007.9 1015.5 1015.4 1054.0 1103.4 1109.1 1143.9
## 2004 1329.0 1477.5 1647.3 1635.6 1529.5 1469.7 1523.4 1563.1 1614.9 1666.7
## 2005 1687.8 1723.8 1773.0 1790.0 1751.7 1938.0 2063.8 2115.8 2133.3 2301.1
## 2006 2677.7 2851.6 2927.3 3610.9 4306.0 3897.7 4170.5 4059.8 4032.9 3998.6
##         Nov    Dec
## 1960  222.8  227.5
## 1961  225.7  226.3
## 1962  230.4  230.5
## 1963  230.6  232.2
## 1964  500.1  453.3
## 1965  523.8  541.4
## 1966  464.4  433.4
## 1967  515.3  551.3
## 1968  458.0  493.4
## 1969  679.4  707.8
## 1970  451.9  435.4
## 1971  406.4  410.7
## 1972  427.8  435.7
## 1973  950.4  959.9
## 1974  608.2  553.2
## 1975  575.1  568.9
## 1976  780.3  767.3
## 1977  650.2  679.1
## 1978  749.3  771.7
## 1979  978.3 1005.6
## 1980  839.5  800.4
## 1981  867.7  869.3
## 1982  884.3  911.4
## 1983  940.4  986.7
## 1984 1085.1 1113.5
## 1985  950.9  962.4
## 1986  915.0  927.6
## 1987 1419.9 1566.6
## 1988 1826.2 1915.3
## 1989 1647.3 1514.1
## 1990 1315.9 1292.7
## 1991 1336.8 1216.4
## 1992 1413.2 1422.5
## 1993 1100.3 1156.4
## 1994 1763.6 1914.4
## 1995 1904.8 1898.7
## 1996 1343.7 1366.9
## 1997 1135.0 1061.5
## 1998  946.9  881.9
## 1999 1065.3 1093.7
## 2000 1258.9 1264.4
## 2001  994.2 1020.8
## 2002 1006.4 1005.7
## 2003 1216.0 1258.5
## 2004 1678.5 1631.1
## 2005 2461.6 2621.7
## 2006 3675.9 3435.8

plot(mcopper)

lambda = BoxCox.lambda(mcopper)
plot(BoxCox(mcopper, lambda))

BoxCox.lambda function suggests using Box-Cox transformation with parameter \(\lambda=0.19\).

Finding a suitable model using auto.arima and trying other models

fit = auto.arima(mcopper, lambda = lambda)
fit

## Series: mcopper 
## ARIMA(0,1,1) 
## Box Cox transformation: lambda= 0.1919 
## 
## Coefficients:
##         ma1
##       0.372
## s.e.  0.039
## 
## sigma^2 estimated as 0.05:  log likelihood=45.05
## AIC=-86.1   AICc=-86.08   BIC=-77.43

Other models which can be tried are ARIMA(2,1,0) and ARIMA(7,1,0):

fit1 = Arima(mcopper, c(2,1,0), lambda = lambda)
fit2 = Arima(mcopper, c(7,1,0), lambda = lambda)

c(fit$aicc, fit1$aicc, fit2$aicc)

## [1] -86.075 -82.512 -80.075

The lowest AICc is for model ARIMA(0,1,1) ( the model which was proposed by auto.arima).

Residuals diagnostics for model ARIMA(0,1,1)

hist(fit$residuals)

tsdisplay(fit$residuals)

Box.test(fit$residuals, lag = 12, type = "Ljung-Box", fitdf = 1)

## 
##  Box-Ljung test
## 
## data:  fit$residuals
## X-squared = 9.48, df = 11, p-value = 0.58

• The residuals look normal at the histogram

• No obvious patterns are visible on the graph of the residuals

• The residuals are white noise

Comparing forecasts of auto.arima and ets

fcast = forecast(fit, h = 20)
plot(fcast)

fit.ets = ets(mcopper)
fcast.ets = forecast(fit.ets, h = 20)
plot(fcast.ets)

Forecasts look differently, but taking into account that the forecasting intervals are very wide, both forecasts provide similar results.

Ex 13

Choose one of the following seasonal time series: hsales, auscafe, qauselec, qcement, qgas.

1. Do the data need transforming? If so, find a suitable transformation.

2. Are the data stationary? If not, find an appropriate differencing which yields stationary data.

3. Identify a couple of ARIMA models that might be useful in describing the time series. Which of your models is the best according to their AIC values?

4. Estimate the parameters of your best model and do diagnostic testing on the residuals. Do the residuals resemble white noise? If not, try to find another ARIMA model which fits better.

5. Forecast the next 24 months of data using your preferred model.

6. Compare the forecasts obtained using ets().

Ex 14

For the same time series you used in the previous exercise, try using a non-seasonal model applied to the seasonally adjusted data obtained from STL. The stlf() function will make the calculations easy (with method="arima"). Compare the forecasts with those obtained in the previous exercise. Which do you think is the best approach?

Ex 15

For your retail time series (Exercise 5 above):

1. develop an appropriate seasonal ARIMA model;

2. compare the forecasts with those you obtained in earlier chapters;

3. Obtain up-to-date retail data from the ABS website (Cat 8501.0, Table 11), and compare your forecasts with the actual numbers. How good were the forecasts from the various models?

Ex 16

1. Produce a time plot of the sheep population of England and Wales from 1867–1939 (data set sheep).

autoplot(sheep)

2. Assume you decide to fit the following model: \[ y_t = y_{t-1} + \phi_1(y_{t-1}-y_{t-2}) + \phi_2(y_{t-2}-y_{t-3}) + \phi_3(y_{t-3}-y_{t-4}) + e_t, \] where \(e_t\) is a white noise series. What sort of ARIMA model is this (i.e., what are \(p\), \(d\), and \(q\))?

ARIMA(3,1,0). \(p=3\), \(d=1\), \(q=0\).

3. By examining the ACF and PACF of the differenced data, explain why this model is appropriate.

sheep %>% diff() %>% ggtsdisplay()

The last significant spike in the PACF is at lag 3, so \(p=3\) is a good choice. Hence an ARIMA(3,1,0) model.

4. The last five values of the series are given below:

Year	1935	1936	1937	1938	1939
Millions of sheep	1648	1665	1627	1791	1797

The estimated parameters are \(\phi_1 = 0.42\), \(\phi_2 = -0.20\), and \(\phi_3 = -0.30\). Without using the forecast function, calculate forecasts for the next three years (1940–1942).

From the equation above \[\begin{align*} \hat{y}_{T+1|T} &= y_{T} + \phi_1(y_{T}-y_{T-1}) + \phi_2(y_{T-1}-y_{T-2}) + \phi_3(y_{T-2}-y_{T-3})\\ &= 1797 + 0.42 (1797-1791) + -0.20 (1791-1627) + -0.30 (1627-1665) \\ &= 1778.00 \\ \hat{y}_{T+2|T} &= \hat{y}_{T+1|T} + \phi_1(\hat{y}_{T+1|T}-y_{T}) + \phi_2(y_{T}-y_{T-1}) + \phi_3(y_{T-1}-y_{T-2})\\ &= 1777.99577 + 0.42 (1777.99577 -1797) + -0.20 (1797-1791) + -0.30 (1791-1627) \\ &= 1718.87 \\ \hat{y}_{T+3|T} &= \hat{y}_{T+2|T} + \phi_1(\hat{y}_{T+2|T}-\hat{y}_{T+1|T}) + \phi_2(\hat{y}_{T+1|T}-y_{T}) + \phi_3(y_{T}-y_{T-1})\\ &= 1718.86855 + 0.42 (1718.86855 - 1777.99577) + -0.20 (1777.99577-1797) + -0.30 (1797-1791) \\ &= 1695.99 \end{align*}\]

5. Now fit the model in R and obtain the forecasts using forecast. How are they different from yours? Why?

forecast(fit, h=3)

##      Point Forecast  Lo 80  Hi 80  Lo 95  Hi 95
## 1940         1778.0 1687.5 1868.5 1639.5 1916.5
## 1941         1718.9 1561.5 1876.2 1478.3 1959.5
## 1942         1696.0 1494.2 1897.8 1387.3 2004.7

They are the same. If you had different values, it is most likely due to rounding errors.

Ex 17

1. Plot the annual bituminous coal production in the United States from 1920 to 1968 (data set ).

autoplot(bicoal)

2. You decide to fit the following model to the series: \[y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \phi_3 y_{t-3} + \phi_4 y_{t-4} + e_t\] where \(y_t\) is the coal production in year \(t\) and \(e_t\) is a white noise series. What sort of ARIMA model is this (i.e., what are \(p\), \(d\), and \(q\))?

ARIMA(4,0,0). \(p=4\), \(d=0\), \(q=0\).

3. Explain why this model was chosen using the ACF and PACF.

ggtsdisplay(bicoal)

The PACF has a signifcant spike at lag 4, but none at higher lags.

4. The last five values of the series are given below.

Year	1964	1965	1966	1967	1968
Millions of tons	467	512	534	552	545

The estimated parameters are \(c = 162.00\), \(\phi_1 = 0.83\), \(\phi_2 = -0.34\), \(\phi_3 = 0.55\), and \(\phi_4 = -0.38\). Without using the forecast function, calculate forecasts for the next three years (1969–1971).

\[\begin{align*} \hat{y}_{T+1|T} &= c + \phi_1 y_{T} + \phi_2 y_{T-1} + \phi_3 y_{T-2} + \phi_4 y_{T-3} \\ &= 162.00 + 0.83 *545 -0.34 *552 + 0.55 *534 -0.38 *512 \\ &= 527.63 \\ \hat{y}_{T+2|T} &= c + \phi_1 \hat{y}_{T+1|T} + \phi_2 y_{T} + \phi_3 y_{T-1} + \phi_4 y_{T-2} \\ &= 162.00 + 0.83* 527.63 -0.34 *545 + 0.55 *552 -0.38* 534 \\ &= 517.19 \\ \hat{y}_{T+3|T} &= c + \phi_1 \hat{y}_{T+2|T} + \phi_2 \hat{y}_{T+1|T} + \phi_3 y_{T} + \phi_4 y_{T-1} \\ &= 162.00 + 0.83* 517.19 -0.34 * 527.63 + 0.55 *545 -0.38* 552 \\ &= 503.81 \\ \end{align*}\]

2. Now fit the model in R and obtain the forecasts from the same model. How are they different from yours? Why?

forecast(fit, h=3)

##      Point Forecast  Lo 80  Hi 80  Lo 95  Hi 95
## 1969         527.63 459.88 595.38 424.02 631.24
## 1970         517.19 429.00 605.38 382.32 652.07
## 1971         503.81 412.48 595.13 364.13 643.48

Again, any differences between these values and those computed “by hand” are probably due to rounding errors.

Ex 18